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ABSTRACT 

We present results from a fully relativistic Monte Carlo simulation of diffusive 
shock acceleration (DSA) in unmodified shocks. The computer code uses a single 
algorithmic sequence to smoothly span the range from nonrelativistic speeds to 
fully relativistic shocks of arbitrary obliquity, providing a powerful consistency 
check. While known results are obtained for nonrelativistic and ultra-relativistic 
parallel shocks, new results are presented for the less explored trans-relativistic 
regime and for oblique, fully relativistic shocks. We find, for a wide trans-relativ- 
istic range extending to shock Lorentz factors > 30, that the particle spectrum 
produced by DSA varies strongly from the canonical f{p) oc spectrum 
known to result in ultra-relativistic shocks. Trans-relativistic shocks may play 
an important role in 7-ray bursts and other sources and most relativistic shocks 
will be highly oblique. 

Subject headings: ISM: cosmic rays — acceleration of particles — relativistic shock 
waves — gamma-ray bursts 



1. Introduction 

Relativistic shocks, and their associ- 
ated energetic particles resulting from dif- 
fusive particle acceleration, have received 
considerable attention in recent years (e.g., 
Ostrowski 1993; Achterberg et al 2001; 
Ellison & Double 2002; Meli & Quenby 
2003a; Niemiec & Ostrowski 2004). This 
is due primarily to their likely presence in 



extreme space phenomena, in particular 7- 
ray bursts (GRBs) (e.g., Piran 1999; Rees 
2000), and their possible role in produc- 
ing ultra-high-energy cosmic rays (UHE- 
CRs) (e.g., Waxman 2000). While most 
work on diffusive shock acceleration (DSA) 
in relativistic shocks has been restricted to 
the ultra-relativistic regime with magnetic 
fields assumed parallel to the shock nor- 
mal, trans-relativistic shocks are certain 
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to be important in some sources and, in 
general, relativistic shocks will have highly 
oblique magnetic fields. Here we consider 
test-particle difi^usive shock acceleration in 
shocks of arbitrary obliquity^ and with ar- 
bitrary Lorentz factors, 70, ranging from 
nonrclativistic shocks to fully relativistic 
ones. 

Diffusive shock acceleration in rela- 
tivistic, oblique shocks was initially ad- 
dressed by Kirk & Heavens (1989) and Bal- 
lard & Heavens (1991) both analytically 
and numerically, and by Ostrowski (1991), 
using numerical techniques. Ostrowski 
(1993) further investigated oblique, rela- 
tivistic shocks using small perturbations 
superimposed on a uniform magnetic field, 
and Bednarz & Ostrowski (1996) inves- 
tigated the corresponding time-scale for 
test-particle acceleration. Bednarz & Os- 
trowski (1998) determined the energy spec- 
tra of accelerated test particles in oblique, 
ultra-relativistic shocks with results lead- 
ing to a limiting energy spectral index of 
as = 2.3±0.1, independent of shock obliq- 
uity. This result had been anticipated by 
Heavens & Drury (1988) for test-particle 
acceleration in ultra-relativistic, parallel 
shocks. More recent work on oblique, rela- 
tivistic shock acceleration has been pre- 
sented by Meh (2002), Meh & Quenby 
(2003b) and Niemiec & Ostrowski (2004). 

Despite the emphasis on ultra-relativ- 
istic shocks in theoretical work, trans-rela- 
tivistic shocks may play an important role 

^Oblique shocks are those where the angle between 
the upstream magnetic field and the shock normal, 
0BO, is greater than 0°. Parallel shocks are those 
with 6bo = 0°. See Fig. 1 for the shock geometry. 
Everywhere in this paper we use the subscript 
(2) to indicate upstream (downstream) quantities. 



in GRBs, both for internal shocks and for 
the expanding fireball shock presumably 
responsible for the afterglow. Even though 
GRB fireballs may initially have Lorentz 
factors 7o > 100, the internal shocks, cred- 
ited by many as responsible for convert- 
ing the bulk kinetic energy of expansion 
into particle internal energy and hence to 
radiation, may be much slower with 70 's 
of a few. As the fireball expands into the 
ambient interstellar medium, regardless of 
its initial 70, it will slow and pass through 
a trans-relativistic phase before becoming 
nonrclativistic. Existing afterglow obser- 
vations span this phase and understand- 
ing Fermi acceleration in trans-relativistic 
shocks is important for interpreting the ob- 
served emission. 

Considering the importance of oblique, 
trans-relativistic shocks, we have extended 
our well-tested Monte Carlo model of DSA 
(e.g., Elhson, Jones, & Reynolds 1990; 
Jones & Ellison 1991; Ellison & Double 
2002) to include shocks of arbitrary obliq- 
uity and speed. We describe the details of 
the model and present results for shocks 
with parameters not previously addressed 
in published work. In the two extreme 
cases where direct comparisons can be 
made - nonrclativistic and ultra-relativis- 
tic shocks - our results are consistent with 
previous work. In future studies we will 
investigate nonlinear effects in oblique rel- 
ativistic shocks where accelerated particles 
modify the shock structure (e.g., Ellison 
& Double 2002). Here we consider only 
particle acceleration in plane, unmodified 
shocks where effects on the shock structure 
from superthermal particles are ignored. 

In contrast to nonrclativistic shocks, 
the details of particle scattering strongly 
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influence the superthermal particle popula- 
tions produced in trans-relativistic and ul- 
tra- relativistic shocks (e.g., Bednarz & Os- 
trowski 1996, 1998). Unfortunately, these 
details are not known with any reliabil- 
ity so the results of all current models of 
DSA depend on the particular scattering 
assumptions made. We use a simple, pa- 
rameterized model of particle diffusion and 
attempt to describe our procedures in suf- 
ficient detail so readers can clearly see how 
the assumptions influence the particle dis- 
tribution functions. While we make no 
claim that our scattering scheme is more 
realistic than other models, some of which 
are far more complex (e.g., Niemiec & Os- 
trowski 2004), we do believe it adequately 
parameterizes the particle transport. Un- 
til a self-consistent theory of wave-particle 
interactions in relativistic shocks is pro- 
duced, parameterization will be necessary. 

Unlike most other models of DSA, 
we follow particles from thermal energies 
through the injection process to superther- 
mal energies. While not necessary in test- 
particle calculations, some description of 
the injection process from thermal ener- 
gies must be used in self-consistent, non- 
linear calculations. The procedures we use 
to model injection in relativistic shocks are 
identical to those we have used with some 
success in nonrelativistic shocks (e.g., Elli- 
son, Mobius, & Paschmann 1990; Jones & 
Ellison 1991; Baring et al. 1997). 

Finally, while our results are applicable 
to ion acceleration in GRBs and the pos- 
sible production of UHECRs, they are not 
directly applicable to the photon emission 
for two basic reasons. The first is that en- 
ergy budget requirements of GRB models 
generally require extremely efficient con- 



version of the bulk kinetic energy of the 
fireball into radiation. This means that the 
particle acceleration process must be effi- 
cient and, therefore, self-consistent, non- 
linear models must be used instead of the 
test-particle ones we present here. 

Second, the radiation seen from GRBs 
is produced by electrons not protons. Nor- 
mal non-relativistic DSA in a proton- 
electron plasma will always put more en- 
ergy into protons than electrons. In rel- 
ativistic shocks, the fraction of energy 
going into electrons is reduced dramati- 
cally. In order to determine the distribu- 
tion of energy between protons and elec- 
trons for GRBs, a two-component acceler- 
ation model must be done in the nonlinear 
regime. This requires additional assump- 
tions for the electron injection and, to ob- 
tain results even remotely consistent with 
GRB properties, requires the modeling of 
lepton dominated plasmas. Preliminary 
work for lepton dominated parallel shocks 
is given in Double (2003). 
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Fig. 1. — Plane shock geometry where the 
shock lies in the y-z plane. 
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2. Monte Carlo Simulation of DSA 

Many of the details of the model we use 
here have been described previously in the 
context of nonrclativistic, oblique shocks 
(e.g., Ellison, Baring, & Jones 1996; Elli- 
son, Jones, & Baring 1999), or rclativistic, 
parallel shocks (Elhson & Double 2002). 
We refer to this previous work and describe 
in greater detail those aspects that are be- 
ing presented for the first time. 

2.1. Jump conditions 

We let Mo be the flow speed of the un- 
shocked (i.e., upstream) plasma as mea- 
sured in a frame at rest with the shock. 
For convenience, we assume the unshocked 
flow is normal to the shock along the x- 
direction so that Uq = u^q. The shock 
Lorentz factor is then 7o = [1 — (mq/c)^]""^/^, 
where c is the speed of light. Down- 
stream, the flow will, in general, have a 
2;-component as well, i.e., U2 = a/m^2 + ""22 
(sec Fig. 1). For nonrclativistic shocks, 
the jump conditions, i.e., the downstream 
density p2, flow speed U2, pressure P2, and 
magnetic fleld B2, in terms of given up- 
stream values, are well-known. For trans- 
relativistic shocks, however, the jump con- 
ditions are not straightforward and can 
only be determined numerically. In all 
of the examples given here, except where 
specifically noted, we assume high sonic 
and Alfven Mach numbers and use the re- 
sults of Double et al. (2004) to calculate 
the jump conditions. 

The primary jump conditions come 
from the equation of continuity which 

gives the compression ratio r = Uxo/ux2 = 
I2P2/ (7oPo)j and the electromagnetic bound- 
ary condition which determines the an- 



gle the downstream magnetic field makes 
with the shock normal, 9b2- Here, 72 = 
[1 — {u2/c)'^]^^^'^ and note that only the 
components of flow speed along the shock 
normal are used in the deflnition of r. 

2.2. Test-particle power law 

For nonrclativistic shock speeds, the 
test-particle power law from DSA is well- 
known and is given by: 

f{p) d^p oc p-'^^P d^p ; (7TP = 3r/(r - 1) , 

(1) 

where p is momentum and f{p) d^p is the 
rotationally averaged, isotropic number 
density of particles in d^p (Krymsky 1977; 
Axford, Leer, & Skadron 1977; Blandford 
& Ostriker 1978; BeU 1978). For high sonic 
and Alfven Mach number, nonrclativistic, 
unmodifled shocks, the compression ratio 
r ~ 4 and the test-particle power-law index 
(7tp, is also approximately 4.^ One of the 
remarkable aspects of test-particle DSA is 
that for Vp ^ Uq, where Vp is the par- 
ticle speed measured in the local plasma 
frame, axp depends only on r regardless of 
the diffusive properties of the plasma and 
eqn. (1) holds regardless of the shock obliq- 
uity (e.g., Blandford & Ostriker 1978). In 
other words, provided that the particles 
are nearly isotropic simultaneously in the 
upstream, shock, and downstream refer- 
ence frames, eqn. (1) holds. 

For rclativistic shocks, -Uq ~ c and par- 
ticle distributions are not isotropic across 
reference frames. This does not change the 
physics of the shock acceleration mecha- 
nism - particles still gain energy by scat- 

^Note that the energy spectral index for fully rel- 
ativistic particles cje, corresponds to a spectral 
index in momentum phase space = ctb + 2. 
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tering back and forth across the shock - 
but it does make obtaining analytic de- 
scriptions of the acceleration process far 
more difficult because the diffusion approx- 
imation can no longer be made. There 
are two important consequences of this. 
First, Monte Carlo computer simulation 
techniques, which do not need to make any 
assumption concerning the isotropy of the 
particle distributions, become the method 
of choice for studying particle accelera- 
tion (e.g.. Kirk & Schneider 1987; Ellison, 
Jones, & Reynolds 1990; Ostrowski 1991; 
Achterberg et al. 2001; Ellison & Double 
2002; Meli & Quenby 2003a), and second, 
due to the mathematical difficulties, there 
are no analytic solutions to compare com- 
puter simulation results to as there are 
with nonrelativistic shocks. 

2.3. Details of the Monte Carlo 
method 

Monte Carlo techniques, while straight- 
forward in principle, are, in fact, diffi- 
cult in practice with many subtle features 
which can produce errors if not imple- 
mented properly. To minimize the possi- 
bility of error, we employ a code which has 
a single algorithmic sequence regardless of 
7o, ^BO) or scattering parameters. This 
allows us to smoothly span the parame- 
ter space from nonrelativistic flow speeds, 
where analytic solutions exist, through the 
trans-relativistic regime where analytic so- 
lutions do not exist, to ultra-relativistic 
shocks, where canonical results again exist 
for parallel shocks. While a smooth tran- 
sition between known results is not a suffi- 
cient condition to guarantee accuracy, it is 
a necessary condition for correct results. 
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Fig. 2. — Schematic representation of 
pitch-angle diffusion. 
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Fig. 3. — Spherical triangle representation 
of a scattering event where the momentum 
changes from p to p' and a particle's pitch 
angle changes from 9^ to 9^. 
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2.3.1. Particle diffusion and convection 

Our implementation of pitch-angle dif- 
fusion is described in detail in Ellison, 

Jones, & Reynolds (1990) and Ellison, Bar- 
ing, & Jones (1996). We assume that scat- 
tering is elastic and isotropic in the lo- 
cal plasma frame and that the scattering 
mean free path A in the local frame is pro- 
portional to the gyroradius Vg. We simu- 
late small-angle scattering by allowing the 
tip of the particle's fluid- frame momen- 
tum vector p to undergo a random walk 
on the surface of a sphere (see Fig. 2). 
If the particle originally had a pitch an- 
gle 6% = arccos{p • B/|p| |B|}, and after a 
time St undergoes a small change in direc- 
tion of magnitude S9, its new pitch angle, 
9^, is related to the old by 



cos^;^ 



cos 9^ cos 59 + sin 9^ sin 69 cos , 

(2) 

where is the azimuthal angle of the mo- 
mentum change dp measured relative to 
the plane defined by the original momen- 
tum p and B (see Fig. 3). After each 
scattering, a new phase angle around the 
magnetic field, 0^ , is determined from the 
old phase angle, 0^ , by 



(f)Q — (f)% + arcsin 



sin0 sin 69 
sin 9^ 



(3) 



where cos 69 is randomly chosen from 
a uniform distribution between 1 and 
cos^^max, and (f) is randomly chosen from 
a uniform distribution between — tt and 
TT, SO that the tip of the momentum vec- 
tor walks randomly over the surface of a 
sphere of radius p = |p| . 

If the time in the local frame required 
to accumulate deflections of the order of 
90° is identified with the coUision time tr — 



X/vp, Ellison, Jones, & Reynolds (1990) 
showed that 



69^ 



(4) 



where 6t is the time between pitch-angle 
scatterings. We take A proportional to 
the gyroradius = pc/ {eB) (e is the 
electronic charge and B is the local uni- 
form magnetic field in Gaussian units), 
i.e., A = rjiiiipTg, where r/mfp is a mea- 
sure of the "strength" of scattering. The 
strong scattering limit, r^mfp = 1, is called 
Bohm diffusion. Now, setting 6t — Tg/Ng, 
where Ng ^ 1 is the number of gyro- 
time segments 6t, dividing a gyro-period 
Tg = 27irg/vp, we have 



S9^eoc = \J 127r/ {r]rafpNg 



(5) 



and the scattering properties of the medium 
are modeled with the two parameters rjmip 
and Ng. 

An important limitation of this scheme 
is the assumption of elastic scattering 
which implies that the scattering centers 
are frozen into the fluid. This eliminates 
both the possibility of second-order Fermi 
acceleration and the transfer of energy 
from accelerated particles to the back- 
ground plasma via the production and 
damping of magnetic turbulence. We also 
neglect any cross-shock electric potential 
that may exist. 

Large iV^'s mean particles make many 
pitch angle scatterings during a gyro- 
pcriod, each with a small angular devia- 
tion. The size of the angular deviation 
depends on r^mfp. A particular value of 
?7mfp means particles will, on average, scat- 
ter through ~ 90° while traversing a dis- 
tance ~ rjraipTg. A large rj^aip implies weak 
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scattering. This implies, of course, that 
magnetic fluctuations with suflicient power 
exist with correlation lengths on the order 
of Lc = '2nrimfprg/Ng. The smaller Ng be- 
comes, the more the scattering resembles 
"large- angle scattering" where the direc- 
tion of p is randomized in a few interac- 
tions with the background magnetic field. 

For nonrelativistic shocks, A^^ has little 
effect on the results and then only when 
Vp ^ Uq. For parallel shocks, regardless 
of 7o, ?7,iifp is unimportant and the value 
of Ny required for "convergence" increases 
as 7g (Ellison & Double 2002). By conver- 
gence we mean that as A^^ is increased, the 
spectral index asymptotically approaches a 
fixed value.^ 

For oblique shocks ^^mfp and Ng are 
both important for the following reasons. 
The size of r^mfp determines the strength of 
cross-field diffusion since on average, every 
time a particle moves ~ A along B, it will 
shift ~ Vg across B. Cross-field diffusion 
is unimportant in parallel shocks but plays 
a critical role in the injection and accel- 
eration processes in oblique shocks since it 
makes it easier for downstream particles to 
rc-cross the shock into the upstream region 
and be further accelerated. In fact, with 
weak cross-field diffusion (i.e., for ?7mfp ^ 
1), downstream particles will not be able to 
re-cross an oblique, ultra-relativistic shock 
when cos^B2 i$ 1/3- With 9^2 determined 
from the jump conditions of Double et al. 
(2004), we see that weak scattering pre- 

^When Ng is less than the convergent value, we 
divide the time between scatterings 5t = Tg/Ng, 
into subdivisions, Stsuh = St/Ns^h- In this case, 
a particle is moved without scattering for TVsub 
steps and the total number of gyro-time segments 
needed for convergence is NgNguh- 



eludes DSA in essentially all shocks with 
7o > 10 unless they arc strictly parallel. 

We inject particles with a thermal dis- 
tribution far upstream and allow them to 
convect into the shock. ^ A particle is trans- 
lated in the following way. During the time 
Stpi between scatterings (measured in the 
local plasma frame) a particle will gyrate 
in the local frame and the local frame will 
convect relative to the shock. For the infi- 
nite, plane shocks we consider, only motion 
in the x-direction is important and for mo- 
tion in the shock frame we have 

Axsk = Ju[Sxpi - Tg sin (cos - 

cos ] + li^^^ipf , (6) 

where 

5xp{ = ps^ipf cos 9b/ (mp7pf) (7) 

is the x-distance the particle moves in the 
local frame, is the component of mo- 
mentum along B (measured in the local 
frame), mp is the mass of a proton, 7pf = 

[1 - {Vp/C)r"\ 7n = [1 - {u,/cf]-"\ 

where Ut is either uq or M2, and (prg {<Prg) 
is the final (initial) phase of the gyrora- 
dius relative to the z'-axis, i.e., in Fig. 4, 

2.3.2. Lorentz transformations and shock 
crossings 

When a particle crosses the shock 
(where the bulk flow and magnetic fleld 
change direction and magnitude) its pitch 
angle and phase change discontinuously. 
However, until a scattering occurs, it 

^In all of the results presented here, the temper- 
ature of the unshocked upstream plasma is low 
enough that the particle ensemble is nonrelativis- 
tic. 
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Fig. 4. — Figure 14 from Ellison, Baring, 
& Jones (1996) showing the geometry at 

the shock where the magnetic field changes 
direction and magnitude. In the notation 
used here 0Bn,i = 6'bo, 0Bn,i+i = 6*32, -Bi = 
i?o, and = i?2. 



moves smoothly and the magnitude of its 
momentum remains constant in a given 
reference frame. The details of the shock 
crossing are given in Ellison, Baring, & 
Jones (1996) and we reproduce Fig. 14 
from that paper here as our Fig. 4. While 
the geometry described in Ellison, Baring, 
& Jones (1996) is identical to that used 
here, there are important differences in 
the techniques we now employ.^ The two 
most important differences are first, that 
our calculation is now fully relativistic and 
applies for all 70 and second, we no longer 
transform to the de Hoffmann- Teller (H- 
T) frame when crossing the shock. Before 
discussing the H-T frame, we describe (fol- 
lowing the arguments in Ellison, Baring, 
& Jones 1996) the calculations required 
in a shock crossing. For convenience of 
discussion we describe a particle crossing 
from upstream to downstream, although 
everything we present applies equally to 
particles crossing in either direction. 

Having determined that a particle has 
crossed the shock, we first transform its lo- 
cal frame momentum (in this case, its mo- 
mentum measured in the upstream frame, 
subscript "0") to the shock frame: 



Pa;,sk 



(7. 



uO 



+ 1 



Py,sk 



(7«o - 1) —Pzo + 

PyO , 



PxO + 

(8) 
(9) 



^Fig. 4 from Ellison, Baring, & Jones (1996) em- 
ploys notation suitable for modified shocks where 
the field and flow change continuously along the x- 
axis. Since here we treat only unmodified shocks, 
the field and flow change only at the shock at 
a; = 0. 
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Pz,sk 



(7 



uO 



1) 1 — ^ ) PxO + 



Un 



(7.0 -1)1^ 

Uq 



luolpirripUzo 



+ 1 



PzQ + 

(10) 



Here, 7„o = [1 — (^o/c)^]""^/^. Given the 
shock frame momentum, we now perform 

a transformation to the downstream frame. 
Note that the values of u^, u^, and 7„ 
change to the downstream values for this 
transformation: 



Px2 



{lu2 - 1) 



ilu2 - 1) 



Ux2 
Ux2Uz2 

ul 



Py2 

Pz2 



lu2lsk'rnpU^2 , 

Py,sk -I 

'Ux2Uz2 



{lu2 - 1) 



ul 



1 Px,sk + 
Pz,sk - (11) 

(12) 

Px,s-k + 



lu2ls]jripUz2 



+ 1 



Pz,sk 



(13) 



where 



and 



7sk = ^[Pt,sk/("^PC)]2 + 1 (14) 



Pt,sk = y^p',sk + pLk + p',sk ■ (15) 

Simultaneously with these transforma- 
tions, the new phase and pitch angle of the 
particle relative to the downstream mag- 
netic field must be computed from the old 
phase and pitch angle measured against 
the upstream magnetic field. This is de- 
scribed in Fig. 4 where, in our current no- 
tation, ©Bn,i = ^BO) ©Bn,j+1 = ^B2, -Bj = 



So, and Bi^i = B2, since we only consider 
unmodified shocks. Referring to Fig. 4, the 
component of the old momentum (i.e., the 
momentum before crossing the shock) in 
the y-direction is 



Py^P±osm(f)% , 



(16) 



and the component along the 2;'-direction 



IS 



Pz' = P±o cos , 



(17) 



where p±o is the component of momentum 
(measured in the upstream frame) perpen- 
dicular to -Bo, and z' is perpendicular to the 
y-Bo plane. The component of momentum 
along the 2;"-direction (i.e., the axis per- 
pendicular to the y-B2 plane) is: 

Pz" = Pz' cos AGbii - Pbo sin AGbii , (18) 

where A©Bn — ^B2 — ^bo- The total mo- 
mentum perpendicular to the new (i.e., 
downstream) magnetic field B2 is 



P±N = \IpI + pI" ■ 



(19) 



Finally, the new phase around B2 is given 

by 



arctan 



Py_ 

Pz" 



(20) 



The two frame transformations and the 
pitch angle and phase transformations 
occur in the instant the particle crosses 
the shock between scatterings. The fact 
that we include in the transformations the 
changing magnetic field strength and di- 
rection and the changing flow speed and 
direction, fully takes into account the u x 
B electric field. No transformation to the 
H-T frame is required. Furthermore, there 
is no need to assume the conservation of 
magnetic moment when crossing the shock. 
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2.3.3. de Hoffmann- Teller frame trans- 
formations 

The de Hoffmann-Teller frame is the 
frame where the u x B electric field is zero. 
For the geometry shown in Fig. 1, the H- 
T frame moves in the negative z-direction 
with speed i'ht = 70^0 tan ^bo- As men- 
tioned above, this frame has been used for 
oblique shocks, particularly with guiding 
center approximations, in attempts to sim- 
plify calculations of the energy change a 
particle experiences in crossing the shock 
(e.g., EUison, Baring, & Jones 1996; Meli 
& Quenby 2003b). Viewed from the shock 
frame, it appears that particles drifting 
in the shock layer experience an energy 
change from the u x B field which is 
not included in normal DSA where par- 
ticles gain energy by scattering between 
the converging upstream and downstream 
plasmas. However, as was evident in the 
original formulation of nonrelativistic DSA 
(e.g., Blandford & Ostriker 1978), where it 
was shown that the power-law index was 
independent of ^^boj shock drift accelera- 
tion must be part of DSA. 

From a particle point of view, particles 
gyrating across an oblique shock can, de- 
pending on their pitch angle, gyrate back 
and forth several times before being con- 
vected downstream or reflected. The en- 
ergy change that occurs in this gyration 
can be viewed as that obtained by the par- 
ticle in the u x B electric field as it moves 
a distance larger than a gyroradius in the 
shock layer. However, since the u x B 
field is zero in the H-T frame and noth- 
ing fundamental changes with a simple 
frame transformation, the energy change 
must be caused by some effect which does 



not depend explicitly on the electric field. ^ 
In fact, as particles gyrate in the shock 
layer they continually cross and re-cross 
the shock and receive repeated energy 
changes from the frame transformations in 
the converging upstream and downstream 
flows. As long as enough scattering occurs 
to maintain near isotropic distributions in 
all frames, this is just the standard DSA 
process and gives exactly the same result 
as obtained from explicitly including the u 
X B electric field. 

For scatter-free propagation, shock- 
drift acceleration can be viewed as an in- 
dependent mechanism. In principle, it is 
possible, depending on the pitch angle a 
particle has in the upstream region, for 
a particle to gain a large amount of en- 
ergy as it gyrates in the shock layer (e.g.. 
Decker 1988). For some pitch angles, up- 
stream particles can even gyrate into the 
shock layer and return upstream, all with- 
out scattering. While in principle it is 
possible for a particular particle to gain a 
large amount of energy in this fashion, the 
range in pitch angles that result in large 
energy gains is extremely small and only 
a small fraction of all particles gain ener- 
gies much in excess of that from a single 
shock compression. Scatter-free propaga- 
tion, even for nonrelativistic shocks, does 
not result in a canonical power law but 
depends strongly on ^bo and the energy of 
the injected particles.^ If elastic scatter- 
ing is included, at any level as long as it 

^We repeat that our work ignores any effects from 
a cross-shock potential which, of course, cannot 
be transformed away. We also ignore any other 
plasma effects which might occur in the shock 
layer such as large amplitude waves or shock surf- 
ing. 

^For thermal particles in high Mach number shocks, 
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is strong enough to drive the distributions 
to isotropy, the process reverts to standard 
DSA. 

In nonrelativistic shocks, a H-T frame 
can be found unless 6*60 — 90° and a H-T 
transformation may be useful in some ap- 
plications. However, a H-T frame with 
speed < c is excluded when tan^Bo > 
c/(mo7o)? making this technique essen- 
tially useless for ultra- relativistic, oblique 
shocks. 

Fortunately, as we have described 
above (and shown in Ellison, Jones, & Bar- 
ing 1999), there is no need to transform to 
the H-T frame and the effects of the u x 
B electric field can be included in the par- 
ticle translation. This allows us to model 
relativistic shocks of arbitrary obliquity. 

2.3.4- Probability of return and distribu- 
tion function 

The particle spectrum produced by 
the Fermi mechanism comes about from 
the average energy gain per shock cross- 
ing combined with the probability that a 
particle will make some number of cross- 
ings. The energy gain per crossing is given 
by the Lorentz transformations discussed 
above. The number of crossings a parti- 
cle makes on average before escaping from 
the shock, is determined by the probability 
that a particle, once downstream, will be 
able to scatter back upstream. This can 
be determined in the Monte Carlo code 
by simply following a particle as it moves 
in the downstream region and assuming. 



no particles will be accelerated beyond simple 

compression by scatter-free propagation in highly 
oblique shocks (e.g., Ellison, Baring, & Jones 
1995). 



once it reaches some number of diffusion 
lengths downstream from the shock that 
it has a vanishing probability of diffusing 
back upstream. A more efficient way of 
doing this is described in Appendix A3 
of Ellison, Baring, & Jones (1996). Once 
a particle becomes isotopic in the down- 
stream frame, the probability of it return- 
ing to some point in the downstream flow 
is 

P^=('-^X. (21) 

This expression is relativistically correct 
(e.g.. Peacock 1981) and applies for any 
shock obliquity. 

For relativistic shocks, the distribution 
function f{p) is calculated in the follow- 
ing way. As downstream particles leave the 
shock system, either from the probability 
of return test or from a set number of diffu- 
sion lengths downstream, they are binned 
in momentum space and f{p) is formed us- 
ing the shock- frame pt,sk- All of our dis- 
tributions with 7o > 2 are calculated this 
way. 

We note that in our previous work (e.g., 
Jones & Ellison 1991; Elhson & Double 
2002) we calculated the distribution at spe- 
cific locations (i.e., specific y — z planes) by 
summing {pt,sk/\Px,sk\) each time a particle 
crossed the plane. This produces an omni- 
directional distribution and is useful for 
simulating what a detector will measure at 
a particular location relative to the shock. 
There is an implicit assumption of isotropy 
with this technique which can be satisfied 
for nonrelativistic shocks and, in nonrel- 
ativistic shocks, the omni-directional dis- 
tribution and the one determined from 
particles leaving downstream are equiva- 
lent as long as the y — z plane is well 
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downstream and all particles are isotropic. 

For relativistic shocks, anisotropics persist 
and the two methods give somewhat differ- 
ent results depending on where the omni- 
directional flux is calculated. 




0.01 0.1 1 10 

P/ClTlpC) 



Fig. 5. — Distribution functions, /(p), ver- 
sus momentum, p, for various input condi- 
tions as labeled. Note that we plot p^f{p) 
in the top panel and p^ '^^f{p) in the bot- 
tom. Spectra with the expected nonrela- 
tivistic, test-particle power-law index will 
be horizontal in this representation. The 
distributions are omni-directional and are 
calculated in the shock frame, downstream 
from the shock. 



3. Results 

3.1. Non-relativistic shocks 

In Fig. 5 we show f{p) versus p for non- 
relativistic shocks (i.e., Uq = 3000 and 10'^ 
km s~^). The spectra are omni-directional 
measured downstream from the shock in 
the shock reference frame. The important 
features in these plots are: 

(i) In the top panel, the distributions 
all become power laws with index a ~ 
(7tp = 3r/(r — 1) ~ 4 when Vp ^ uq. 
This conflrms that our technique, which 
does not make a transformation to the H- 
T frame, properly accounts for the energy 
gain in the shock layer without explicitly 
including the u x B electric field. To em- 
phasize this point, we show in the bottom 
panel the spectrum from a shock with a low 
Alfven Mach number Ma — 4.6. The pa- 
rameters for this shock are: proton density 
= 1 cm~^, unshocked proton temperature 
= 10^ K, unshocked electron temperature 
= 10^ K,s and Bq = 3 x 10"^ G. In this 
case, the magnetic field is strong enough 
to produce a significant change in down- 
stream flow direction (i.e., 9^2 — 10.1°) 
and to modify the jump conditions giving 
a compression ratio of r ~ 3.43. Here, we 
plot p^-^^fip), where 4.23 = 3r/(r - 1), 
and the horizontal power law confirms that 
we get the correct cttp. The dynami- 
cally important magnetic field makes this 
a stronger test of our translation technique 
than the cases in the top panel where Bq = 
10-6 G. For the ^bo = 60° shocks in the top 
panel, the sonic Mach number Ms — 81, 
Ma 4600 and 9^2 ^ 1-4 x 10"^ degrees. 

'The electron temperature is used only for deter- 
mining the Mach numbers. 
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(ii) The normalization of the power law 
relative to the thermal peak drops as ei- 
ther r^mfp or ^^Bo is increased. This drop 
occurs before particles become isotropic 
and comes about because particles with 
Vp ^ Uq enter the downstream region di- 
rected predominantly along the shock nor- 
mal. This makes it harder for them to 
scatter back upstream and the difficulty in- 
creases when either r^mfp or ^bo is large. 

Particles gyrate around and along the 
magnetic field. In the absence of cross-field 
diffusion, the greater 6^2, the further par- 
ticles have to move along the field to move 
a given distance normal to the shock. This 
increases the hkelihood that they will con- 
vect downstream without further acceler- 
ation. When ?7infp is large, cross-field dif- 
fusion is less important compared to scat- 
tering along the magnetic field and parti- 
cles again are less likely to scatter back up- 
stream. 

Note that the probability of return 
equation (21), which does not include 
9b2, is still valid even though particles are 
swept downstream more quickly in oblique 
shocks than in parallel ones. As a down- 
stream particle gyrates around the mag- 
netic field, it will cross a particular point 
in the downstream flow more often if the 
field is oblique than if it is parallel. If this 
point is used to calculate Pret, the par- 
ticle will have some probability of escap- 
ing every time it gyrates across this point 
and is more likely to escape downstream 
than would be the case in a parallel shock. 
Nevertheless, the test-particle power law 
index iTtp is obtained once particles be- 
come isotropic since, on average, a particle 
gains more energy per shock crossing in 
an oblique shock than in a parallel one 



because it can gyrate across the shock sev- 
eral times in a single shock crossing event. 
These two effects, more energy per crossing 
but smaller probability of making cross- 
ings, each depend in the same way on 9-B2 
or r^mfp for isotropic particles and combine 
to give the canonical cttp once Vp ^ Uq. 

(iii) The minimum momentum where 
the power-law tail obtains increases with 
increasing 9bo- This is particularly notice- 
able with the ^Bo — 89° example in the top 
panel of Fig. 5 and comes about because, 
to obtain (Ttp, particles must have a speed 
large compared to the effective ffow speed. 
In oblique shocks, this speed is essentially 
the H-T velocity and the momentum where 
(7tp is established increases with ^bo- 

The non-power-law tail coming off the 
thermal distribution in oblique shocks may 
be important in the heliosphere where 
highly oblique interplanetary traveling 
shocks accelerate the thermal solar wind 
(e.g.. Baring et al. 1997). The energetic 
particles observed by spacecraft may not 
show the canonical Fermi power law even 
though DSA is the physical mechanism 
producing the energetic population. 

3.2. Trans-relativistic shocks 

As soon as the flow speed becomes com- 
parable to the speed of light, the power-law 
index starts to depend strongly on ?7mfp and 
^Bo- In Fig- 6 we show results for a trans- 
relativistic shock with — 2 for various 
?7infp and ^bo- These shocks all have r ~ 
3.32 as determined by the jump conditions 
given in Double et al. (2004) and we plot 
P^-^V(p), where 4.29 ~ 3r/(r - 1). The 
spectra are measured in the shock frame 
from escaping particles. In all cases ex- 
cept for the light-weight solid curve in the 



13 




J I I U\ I I 

1 2 3 4 5 
LoQio [p/(mpC)] 



Fig. 6. — Distribution functions, f{p), ver- 
sus momentum, p, for various input con- 
ditions as labeled. Note that we plot 
where 4.29 ~ 3r/(r - 1) = cttp, 
i.e., the spectral index expected for nonrel- 
ativistic shocks with r = 3.32. The light- 
weight solid curve in the top panel uses a 
small value of Ng. These are trans- relativ- 
istic shocks with 70 = 2 and show spectra 
both flatter and steeper than cttp- 



top panel, the value of Ng is large enough 
to produce convergent results. The heavy- 
weight solid curve in the top panel is the 
parallel shock case and has a ~ 4.1. The 
other curves in the top panel show results 
with ^Bo = 60° for various r^mfp- A power 
law is obtained in each case, but there is a 
strong steepening of the power-law portion 
of the spectrum with increasing ?7mfp- 

In the bottom panel of Fig. 6 we show 
the highly oblique case ^bo = 89° for vari- 
ous ?7mfp- The spectra are similar to those 
with ^Bo = 60° except the steepening with 
increasing ^^mfp is greater. 

As 7o increases, the influence of Ng in- 
creases, but it is still small for 70 = 2. The 
light-weight solid curve in the top panel of 
Fig. 6 was calculated with the same pa- 
rameters as the dotted curve except that 
Ng ^ 20 rather than 1200 for the dot- 
ted curve. The small Ng produces a slight 
steepening of the power law. 

3.3. Fully relativistic shocks 

In Fig. 7 we show results with 70 = 10 
for various r/mfp and ^bo- These shocks all 
have r ~ 3.02 (cttp — 4.49) but we have 
plotted where 4.23 is the index 

known to result for fully relativistic, par- 
allel shocks with strong scattering (e.g., 
Ballard & Heavens 1991; Bcdnarz & Os- 
trowski 1998; Achterberg et al. 2001; Elli- 
son & Double 2002). In all cases, the value 
of Ng is large enough to produce conver- 
gent results. The solid curve in the top 
panel confirms that we obtain the canoni- 
cal result for ^^bo = 0°. 

For 9bo > 0°, the results show the same 
general trend as those for the trans-rela- 
tivistic shocks in Fig. 6. The top panels of 
Figs. 6 and 7 are similar except the steep- 
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Fig. 7. — Distribution functions, f{p), ver- 
sus momentum, p, for shocks with 70 = 10. 
Note that the canonical power law index 
for ultra-relativistic shocks is cr ~ 4.23 and 
we plot to emphasize this result. 



10 



1 r 



0.1 



10 



0.1 



70=20 

Solid: ')„,p=l. 630=0° 
Dashed: 7,^,^=1, 660=60° 
Dotted: 7,^,p=3, 9^0=60° 
Dot-dashed: 7)^fp=3, 9gQ=89° 




2 3 4 5 
Log^o p/(m c) 



Fig. 8. — Distribution functions, f{p), ver- 
sus momentum, p, for shocks with 70 = 20 
(top panel) and 70 = 30 (bottom panel). 
We plot p^''^^f{p) to emphasize the ^bo = 
0° result. 
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ening for 70 = 10 is a less dramatic func- 
tion of ?7mfp- The same is true for the bot- 
tom panels where 9bo = 89°. 
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Fig. 9. — Distribution functions, f{p), ver- 
sus momentum, p, for shocks with 70 — 20. 
The parameter A''^ is varied while all other 
input parameters are held constant. 

The computation time required to pro- 
duce convergent results (i.e., results with 
large enough Ng) increases with 70 so we 
show a limited number of examples with 
7o > 10. Fig. 8 shows results for 70 = 20 
(top panel) and 70 = 30 (bottom panel) 
and these examples, combined with our 
7o = 10 result, show a clear trend. First 
of all, the standard a ~ 4.23 result for 
^Bo = 0° is obtained in all cases. For 
^Bo > 0°, we get results very similar to 
those for 70 = 10 with a general reduction 
in the variation caused by varying 9bo or 
?7mfp- It is likely that this trend will con- 
tinue, supporting the assertion that ultra- 
relativistic shocks produce power laws with 



a ~ 4.23 independent of ^^bo or ?7mfp- This 
assumes, of course, that Ng is large enough 
to produce convergence. 

The effect of varying Ng when 70 = 20 
and ^Bo = 0° is shown in Fig. 9 (simi- 
lar results are obtained for 70 = 10 and 
30). This parameter now dramatically in- 
fluences the spectrum. The step-function 
effect in the spectrum has been noted by 
several authors in the context of large- 
angle scattering (e.g., Quenby & Lieu 1989; 
Ellison, Jones, & Reynolds 1990; Baring 
1999) (see also Vietri 1995). Fig. 9 shows 
that this effect emerges smoothly from the 
canonical cr = 4.23 power law as Ng is low- 
ered and scattering becomes coarser. No- 
tice that the momentum of the first peak 
above thermal energies depends on Ng. For 
Ng = 100 this peak occurs at ~ lO^m^c 
rather than •jQiripC because the curves in 
Fig. 9 (and all others in this paper) are 
calculated in the shock frame. Plotted in 
the downstream plasma frame, this peak is 
at ~ ^QiripC as expected. 

In Fig. 10 we show the effect of coarse 
scattering (small Ng — 40) as a function 
of 77infp for 7o = 20 and ^bo = 60°. The 
step-like structure remains as the spectra 
steepen with increasing ?7mfp, although the 
steps begin to smooth some. 

3.4. Pitch-angle distributions 

The anisotropic nature of particles in 
relativistic shocks produces the strong de- 
pendencies we've seen on ^boj 'Hm.tp-i ^^^d 
Ng. The main difference between nonrel- 
ativistic and relativistic shocks is shown 
in Fig. 11, where we compare the pitch- 
angle distributions of shock crossing parti- 
cles. The curves are calculated by sum- 
ming the quantity Px,sk/pt,sk as particles 
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Fig. 10. — Distribution functions with 
7o = 20 for small A^^^ = 200 and various 
6*60 and r^mfp, as indicated. 

cross the shock and only supcrthcrmal par- 
ticles are included in the plots. These are 
shock-frame values of momentum and the 
dashed curve from our nonrelativistic ex- 
ample (Fig. 5) shows the flux weighting 
characteristic of an isotropic distribution, 
i.e., particles cross the shock with a fre- 
quency proportional to their component of 
velocity normal to the shock. 

The relativistic shocks show highly 
anisotropic pitch-angle distributions with 
most particles crossing into the down- 
stream region at an oblique angle. For 
a given set of parameters 70, ^bo, ^mfp, and 
Ng, a particular pitch-angle distribution 
results and this determines the average 
momentum gain in a shock crossing. This 
combined with the probability to make 
some number of crossings determines the 
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Fig. 11. — Pitch-angle distributions for 
shock crossing particles in shocks ranging 
from nonrelativistic to fully relativistic, as 

indicated. A positive Px,sk indicates motion 
toward the downstream direction. The 
area under each curve is normalized to one. 
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spectrum. Unlike for nonrelativistic shocks 
with their isotopic distributions, no simple 
analytic expression has been found for a 
once flows become trans-relativistic. 
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Fig. 12. — Pitch-angle distributions for the 
shocks shown in Fig. 7. The area under 
each curve is normalized to one. 

In Fig. 12 we show the pitch-angle dis- 
tributions for the examples shown in Fig. 7. 

It's clear that relatively subtle changes in 
the pitch-angle distribution causes fairly 
large changes in /(p). The bottom panel 
for ^Bo = 89° shows the trend as r/mfp in- 
creases; the larger r^mfp, the more peaked 
the distribution is at Px,sk/Pt,sk ~ 0.6 with 
a smaller fraction of particles crossing back 
into the upstream region. 
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Fig. 13. — Pitch-angle distributions for 
shocks with parameters similar to those in 
Fig. 9. The area under each curve is nor- 
malized to one. 

Finally, in Fig. 13 we show pitch-angle 
distributions as Ng is varied as in Fig. 9. 
When Ng is small, the pitch-angle distribu- 
tion is strongly peaked in the forward di- 
rection and this is responsible for the step- 
hke nature of /(p). As Ng is increased, the 
distribution smoothly moves to a configu- 
ration similar to those shown in Fig. 12. 

We illustrate this in Fig. 14. When 
fp ^ Mo, flux weighting causes most parti- 
cles crossing into the upstream region from 
downstream to be directed close to the 
shock normal, i.e., with pitch angles near 
180°. A large Ng means a particle will in- 
teract after a short bt with a small 5^max 
and is less likely to broaden its pitch an- 
gle much before being swept back across 
the shock. This is illustrated in Fig. 14 
with the particle trajectory labeled (A). If 
Ng is small, (^^max is larger, the particle 
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Fig. 14. — Schematic representation of 
particles crossing a shock with large and 
small values of A^^^. 

changes its pitch angle by a larger amount 
on average, and there is a greater chance 
the particle will cross back into the down- 
stream region with a flatter trajectory as 
shown with the trajectory labeled (B). On 
average, particles crossing as in (A) with 
a large Ng receive less energy per crossing 
than those crossing as in (B), where Ng is 
small. 

4. Conclusions 

We have presented a comprehensive 
study of test-particle, diffusive shock accel- 
eration in plane shocks of varying Lorentz 
factor 7o, obliquity 9^0, and scattering 
properties parameterized by r^^fp and Ng. 
We use a Monte Carlo computer simula- 
tion which smoothly spans the parameter 
space from nonrelativistic shocks through 
trans-relativistic ones to fully relativistic 
shocks, without requiring a transforma- 
tion to the de Hoffmann- Teller frame or 
making any assumptions concerning the 
magnetic moment of shock crossing parti- 



cles. Due to the Lorentz transformation of 
the magnetic field, relativistic shocks will 
have superluminal H-T speeds for virtually 
all obliquities, making transformations to 
the H-T frame of limited use. 

For nonrelativistic shock speeds, the 
test-particle Fermi power law has a well- 
known analytic form which is used to test 
our computer simulations. For trans-rela- 
tivistic and fully relativistic shocks, how- 
ever, analytic results are rare or non- 
existent although several independent cal- 
culations provide a canonical power-law for 
ultra-relativistic particles in parallel, ultra- 
relativistic shocks, i.e., /(p) oc p-4-23±o.i_ 
We obtain this result and the fact that 
our simulation has a single algorithmic se- 
quence spanning all 70's and 6'bo's provides 
some measure of assurance that our results 
are accurate in trans-relativistic regions 
where analytical results do not exist. 

In contrast to nonrelativistic shocks, 
DSA in relativistic shocks depends criti- 
cally on the details of how particles inter- 
act with the background magnetic turbu- 
lence. In real shocks, this turbulence must 
be self generated by the shock accelerated 
particles. However, the scattering inter- 
actions responsible for this generation are 
not known with any reliability and particle 
diffusion must be parameterized. Despite 
the complexity and unknown nature of the 
wave-particle interactions, we believe that 
much of the essential underlying physics 
can be modeled by first assuming that the 
scattering mean free path A is proportional 
to the gyroradius Vg and then defining two 
parameters which control the finer scatter- 
ing details. These parameters are ?7mfp in 
the relation A = rjmfpfg, which character- 
izes the "strength" of cross-field scattering. 
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and the "fineness" of pitch-angle scattering 
A''^, where particles pitch-angle scatter af- 
ter a fraction of a gyro-period 5t = Tg/Ng 
within a maximum angle given by eqn. (4). 
Spectra converge to a particular form as Ng 
increases and the convergent value of Ng is 
proportional to 

Our results show that particle distri- 
butions depend strongly on ?7infp and Obo, 
as well as on Ng for relativistic shocks 
when Ng is below the convergent value. 
Most importantly, we show that f{p) de- 
parts significantly from f{p) oc p^'^-'^^ in a 
wide trans-relativistic regime extending to 
above 70 = 30. 

In our examples with 70 = 2 we find 
that the power law is fiatter than cttp = 
3r/(r - 1) (Fig. 6) for ^bo = 0°. That is, 
the power-law index a ^ 4.12 versus o"tp — 
4.29. As long as 9bo = 0°, the results are 
independent of ?7mfp- For 6*30 > 0°, the 
power law steepens as either ^bo or rjmfp is 
increased. For our most extreme 70 = 2 
case, ^Bo — 89° and rj^fy — 6, yielding a ~ 
6.5. 

For our more strongly relativistic ex- 
amples (70 > 10), with convergence val- 
ues of Ng, we obtain the canonical power- 
law index a ~ 4.23 for ^bo = 0°. As with 
7o = 2, the power law steepens when either 
6'bo or r/mfp is increased (Figs. 7 and 8), but 
the steepening becomes a weaker function 
of both ^Bo and r^mfp as 70 increases. Our 
results are consistent with cr 4.23 inde- 
pendently of ^Bo and ?7infp in the ultra-rela- 
tivistic limit as widely reported (e.g., Bed- 
narz & Ostrowski 1998; Achterberg et al. 
2001). We do show, however, that a wide 
range of Lorentz factors exists extending 
up to 7o ~ 30 where a differs significantly 
fi-om 4.23. 



In contrast to the effects of ^bo and 
Vmfp, which diminish with increasing 70, 
the influence of Ng remains as 70 increases. 
Figs. 9 and 10 show that small Ng (large- 
angle scattering) produces large features in 
the spectra as individual shock crossings 
continue to show in f{p) to high energies. 
There is no reason to believe this effect will 
diminish with increasing 70. 

It is not at all obvious what realistic 
values of Ng are or if they are large enough 
to produce convergence. Large values of 
A''^ imply that signiflcant power in mag- 
netic turbulence exists at extremely small 
length scales and it is possible that ac- 
tual plasmas have some lower length-scale 
limit which may be larger than needed 
for convergence. In this case, the spec- 
trum will depend on the effective Ng and 
may be highly variable depending on the 
particle shock parameters. Furthermore, 
the highest momentum particles require 
magnetic turbulence with extremely long 
length scales for resonance. If this turbu- 
lence does not exist, the spectrum will turn 
over at some characteristic momentum de- 
termined by the longest magnetic length 
scale (see Niemiec & Ostrowski 2004, for a 
discussion of such effects). 

Some of the most exotic and interest- 
ing astrophysical objects are likely to har- 
bor relativistic shocks. While it is probable 
that these shocks accelerate particles, the 
spectrum, even in unmodified, ultra-rela- 
tivistic shocks, depends on the unknown 
details of the wave-particle interactions. 
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